program mynormal12, sortpreserve
	version 8.1
	if replay() {
		if (`"`e(cmd)'"' != "mynormal12") error 301
		Replay `0'
	}
	else	Estimate `0'
end

program Estimate, eclass
	syntax varlist [if] [in]		///
		[fweight pweight] [,		///
		noLOg				/// -ml model- options
		noCONStant			///
		Robust				///
		CLuster(varname)		///
		HETero(varlist)			///
		noLRTEST			///
		OFFset(varname numeric)		///
		EXPosure(varname numeric)	///
		svy				///
		Level(integer `c(level)')	/// -Replay- options
		EForm				///
		*				/// -mlopts/svyopts- options
	]

	// mark the estimation sample
	marksample touse

	// check syntax
	if "`svy'" != "" {
		svymarkout `touse'
		local wvar `r(weight)'		// remember weight var

							// NEW block
		svyopts modopts diopts options ,	///
			`options'			///
			`eform'				///
			level(`level')
		if "`s(subpop)'" != "" {
			local subpop `s(subpop)'
			local subopt subpop(`subpop')
		}
	}
	else	local diopts `eform' level(`level')	// <- NEW
	mlopts mlopts, `options'
	local cns `s(constraints)'
	gettoken lhs rhs : varlist
	if "`cluster'" != "" {
		local clopt cluster(`cluster')
	}
	if "`weight'" != "" {
		tempvar wvar
		quietly gen double `wvar' `exp'
		local wgt "[`weight'=`wvar']"
	}
	if "`wvar'" != "" {
		// for initial value calculations
		local awgt "[aw=`wvar']"
	}
	if `"`robust'`cluster'`cns'"' != ""	///
	| "`weight'" == "pweight" {
		local lrtest nolrtest
	}
	if "`log'" != "" {
		local qui quietly
	}
	if "`offset'" != "" {
		local offopt "offset(`offset')"
	}
	if "`exposure'" != "" {
		local expopt "exposure(`exposure')"
	}

	// markout missing values from the estimation sample
	markout `touse' `wvar' `offset' `exposure' `subpop'	// <- CHANGED
	markout `touse' `cluster', strok
	if "`subpop'" != "" {					// NEW block
		// restrict initial value calculations to
		// observations within the subpopulation
		local subtouse `touse' & `subpop' != 0
	}
	else	local subtouse `touse'

	if "`constant'" == "" {
		// initial values
		quietly sum `lhs' `awgt' if `subtouse'		// <- CHANGED
		local mean = r(mean)
		local lnsd = ln(r(sd))+ln((r(N)-1)/r(N))/2
		local initopt init(/mu=`mean' /lnsigma=`lnsd') search(off)

		`qui' di as txt _n "Fitting constant-only model:"
		ml model lf mynormal_lf				///
			(mu: `lhs' =, `offopt' `expopt' )	///
			(lnsigma: `hetero')			///
			`wgt' if `touse',			///
			`log'					///
			`mlopts'				///
			`initopt'				///
			`svy'					///
			`subopt'				/// <- NEW
			nocnsnotes				///
			missing					///
			maximize
		if "`lrtest'" == "" {
			local contin continue search(off)
		}
		else {
			tempname b0
			mat `b0' = e(b)
			local contin init(`b0') search(off)
		}
		`qui' di as txt _n "Fitting full model:"
	}

	// fit the full model
	ml model d2 mynormal_d2				///
		(mu: `lhs' = `rhs',			///
			`constant' `offopt' `expopt'	///
		)					///
		(lnsigma: `hetero')			///
		`wgt' if `touse',			///
		`log'					///
		`mlopts'				///
		`robust'				///
		`clopt'					///
		`contin'				///
		`svy'					///
		`modopts'				/// <- NEW
		missing					///
		maximize

	if "`hetero'" == "" {
		ereturn scalar k_aux = 1
	}
	else	ereturn scalar k_aux = 0
	ereturn local cmd mynormal12

	Replay , `diopts'				// <- CHANGED
end

program Replay
	syntax [, Level(integer `c(level)') * ]
	// use svyopts to get display options
	svyopts noopts diopts, `options'		// <- NEW
	// anything in `noopts' will cause an error
	mlopts noopts, `noopts'				// <- NEW
	if `e(k_aux)' {
		local sigma diparm(lnsigma, exp label("sigma"))
	}
	ml display , level(`level') `sigma' `diopts'	// <- CHANGED
end
